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ABSTRACT 

This paper introduces the multiband periodogram , a general extension of the 
well-known Lomb-Scargle approach for detecting periodic signals in time-domain 
data. In addition to advantages of the Lomb-Scargle method such as treatment 
of non-uniform sampling and heteroscedastic errors, the multiband periodogram 
significantly improves period finding for randomly sampled multiband light curves 
(e.g., Pan-STARRS, DES and LSST). The light curves in each band are modeled 
as arbitrary truncated Fourier series, with the period and phase shared across all 
bands. The key aspect is the use of Tikhonov regularization which drives most 
of the variability into the so-called base model common to all bands, while fits 
for individual bands describe residuals relative to the base model and typically 
require lower-order Fourier series. This decrease in the effective model complexity 
is the main reason for improved performance. We use simulated light curves and 
randomly subsamplcd SDSS Stripe 82 data to demonstrate the superiority of 
this method compared to other methods from the literature, and find that this 
method will be able to efficiently determine the correct period in the majority 
of LSST’s bright RR Lyrae stars with as little as six months of LSST data. A 
Python implementation of this method, along with code to fully reproduce the 
results reported here, is available on GitHub. 

Subject headings: methods: data analysis — methods: statistical 


1. Introduction 


Many types of variable stars show periodic flux variability (Eyer & Mowlavi 2008). 


Periodic variable stars are important both for testing models of stellar evolution and for 
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using such stars as distance indicators (e.g., Cepheids and RR Lyrae stars). One of the first 
and main goals of the analysis is to detect variability and to estimate the period and its 
uncertainty. A number of parametric and non-parametric methods have been proposed to 
estimate the period of an astronomical time series (e.g., Graham et al. 2013, and references 
therein). 

The most popular non-parametric method is the phase dispersion minimization (PDM) 


introduced by Stellingwerf (1978). Dispersion per bin is computed for binned phased light 


curves evaluated for a grid of trial periods. The best period minimizes the dispersion per bin. 
A similar and related non-parametric method that has been recently gaining popularity is the 
Supersmoother routine (Reimann 1994). It uses a running mean or running linear regression 
on the data to fit the observations as a function of phase to a range of periods. The best 
period minimizes a hgure-of-merit, adopted as weighted sum of absolute residuals around 
the running mean. Neither the Supersmoother algorithm nor the PDM method require a 
priori knowledge of the light curve shape. We note that neither method produces posterior 
probability distribution for the period but rather a single point estimate. 

The most popular parametric method is the Lomb-Scargle periodogram, which is dis¬ 
cussed in detail in Section [2} The Lomb-Scargle periodogram is related to the y 2 for a 
least-square fit of a single sinusoid to data and can treat non-uniformly sampled time series 
with heteroscedastic measurement uncertainties. The underlying model of the LombScargle 
periodogram is nonlinear in frequency and basis functions at different frequencies are not 
orthogonal. As a result, the periodogram has many local maxima and thus in practice the 


global maximum of the periodogram is found by grid search (for details see, e.g. Ivezic et al. 


2014). A more general parametric method based on the use of continuous-time autoregressive 


moving average (CARMA) model was recently introduced by Kelly et al. (2014). CARMA 
models can also treat non-uniformly sampled time series with heteroscedastic measurement 
uncertainties, and can handle complex variability patterns. 

A weakness of the above methods is that they require homogeneous measurements - for 
astronomy data, this means that successive measurements must be taken through a single 
photometric bandpass (filter). This has not been a major problem for past surveys because 


measurements are generally taken through a single photometric filter (e.g. LINEAR, Sesar 


et al. 2011), or nearly-simultaneously in all bands at each observation (e.g. SDSS, Sesar 


et al. 2010). For the case of simultaneously taken multiband measurements, Siiveges et al. 


(2012) utilized the principal component method to optimally extract the best period. Their 


method is essentially a multiband generalization of the well-known two-band Welch-Stetson 


variability index (Stetson 1996). Unfortunately, when data in each band are taken at different 


times, such an approach in not applicable. In such cases, past studies have generally relied 
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on ad hoc methods such as a majority vote among multiple single-band estimates of the 


periodogram (e.g., Olnseyi et ah 2012). 


For surveys that obtain multiband data one band at a time, such as Pan-STARRS 


(Kaiser et ah 2010) and DES (Flaugher 2008), and for future multicolor surveys such as 


LSST (Ivezic et ah 2008), this ad hoc approach is not optimal. In order to take advantage of 
the full information content in available data, it would be desirable to have a single estimate 
of the periodogram which accounts for all observed data in a manner which does not depend 
on the underlying spectrum of the object. We propose such a method in this paper. 

The proposed method is essentially a generalization of the Lomb-Scargle method to 
multiband case. The light curves in each band are modeled as arbitrary truncated Fourier 
series, with the period, and optionally the phase, shared across all bands. The key aspect en¬ 


abling this approach is the use of Tikhonov regularization (discussed in detail in Section 4.3) 
which drives most of the variability into the so-called base model common to all bands, while 
fits for individual bands describe residuals relative to the base model and typically require 
lower-order Fourier series. This decrease in effective model complexity is the main reason 
for improved performance. 

The remainder of the paper is organized as follows: in Section [2] we provide a brief 
review of least-squares periodic fitting, and in Section [3] derive the matrix-based formalism 
for single-band periodic analysis used through the rest of this work. Section [4] introduces 
several extensions and generalizations to the single-band model that the matrix formalism 
makes possible, including floating mean models, truncated Fourier models, and regularized 
models. In Section [5j we use the ideas behind these extensions to motivate the multiband 
periodogram , and show some examples of its use on simulated data. In Section [6] we apply this 
method to measurements of 483 RR Lyrae stars first explored by Sesar et al. (2010, hereafter 
S10), and in Section[7]explore the performance of the method for simulated observations from 
the LSST survey. We conclude in Section [8} 


2. Brief Overview of Periodic Analysis 


The detection and quantification of periodicity in time-varying signals is an important 
area of data analysis within modern time-domain astronomical surveys. For evenly-spaced 
data, the periodogram , a term coined by Schuster (1898), gives a quantitative measure of the 
periodicity of data as a function of the angular frequency c o. For data {yk}k=i measured at 
equal intervals t k = t 0 + kAt, Schuster’s periodogram, which measures the spectral power as 
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a function of the angular frequency, is given by 

N 


CM = 4 


X] y* e 


icjtk 


k =1 


(i) 


and can be computed very efficiently using the Fast Fourier Transform. 


Because astronomical observing cadences are rarely so uniform, many have looked at 
extending the ideas behind the periodogram to work with unevenly-sampled data. Most 
famously, Lomb (jl976) and Scargle (1982) extended earlier work to define the normalized 
periodogram-. 


Pn(w) — 


2 K 


E k(Vk - y ) COS u(t k - r)] 2 E k {Vk - y ) sinuj(t k - r)] 2 


Efc COS 2 U)(t k - T ) 


E fc sin 2 o;(t fc - r) 


( 2 ) 


where y is the mean and V y is the variance of the data {y k }, and r is the time-offset which 
makes P/v(ca) independent of a translation in t (see "Press et al. 2007, for an in-depth discus¬ 
sion). Lomb (1976) showed that this time-offset has a deeper effect: namely, it makes Pn 


identical to the estimate of harmonic content given a least-squares fit to a single-component 
sinusoidal model, 

d(t) = j4sin(u;£ + 0). (3) 

This long-recognized connection between spectral power and least squares fitting methods 
was solidified by Jaynes (1987|), who demonstrated that the normalized periodogram of Lomb 
and Scargle is a sufficient statistic for inferences about a stationary-frequency signal in the 


presence of Gaussian noise. Building on this result, Bretthorst (1988) explored the extension 


of these methods to more complicated models with multiple frequency terms, non-stationary 
frequencies, and other more sophisticated models within a Bayesian framework. 


While the important features of least squares frequency estimation via Lomb-Scarglc 
periodograms have been discussed elsewhere, we will present a brief introduction to the 
subject in the following section. The matrix-based formalism we develop here will make 
clear how the method can be extended to more sophisticated models, including the multiband 
periodogram proposed in this work. 


3. Standard Least Squares Spectral Fitting 

In this section we present a brief quantitative introduction to the least squares fitting 
formulation of the normalized periodogram of Equation ([2]). We denote N observed data 
points as 


D {^ki yki < Pc}fc=i 


( 4 ) 
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where t k is the time of observation, y k is the observed value (typically a magnitude), and a k 
describes the Gaussian errors on each value. Without loss of generality we will assume that 
the data yk are centered such that the measurements within each band satisfy 

WkVk 
k W k 

where the weights are w k = <J k 2 - 



3.1. Stationary Sinusoid Model 


The normalized periodogram of Equation (J2]) can be derived from the normalized y 2 of 
the best-fit single-term stationary sinusoidal model given in Equation ([3]) . To make the prob¬ 
lem linear, we can re-express the model in terms of the parameter vector 6 = [A cos (p, A sin cp] 
so that our model is 

y(t\ui, 6) = 9i sin(u;£) + 9 2 cos(cut). (6) 

For a given u, the maximum likelihood estimate of the parameters 9 can be found by mini¬ 
mizing the y 2 of the model, which is given by 

[Vk - y(t k \u,9)} 2 


x 2 M = X] 


cn 


( 7 ) 


For the single-term Fourier model, it can be shown (see, e.g. Ivezic et ah 2014) that 

XmmM = Xot 1 - PnH] 


( 8 ) 


where Pn(uj) is the normalized periodogram given in Equation ([2j and yg is the reference 
y 2 for a constant model, which due to the assumption in Equation ([5]) is simply Xo = 
Efc(Ww-) 2 - 


3.2. Matrix Formalism 

The expressions related to the stationary sinusoid model can be expressed more com¬ 
pactly by defining the following matrices: 



sin(o;£i) 

cos(a;ti) 
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0 • 

• 0 

AG = 

sin^G) 

cos(Lot 2 ) 
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sin (a An) 

COS (0Jt]y) 
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(9) 
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With these definitions, the model in Equation (|6]) can be expressed as a simple linear product, 
y(t\u, 6) = X u 9, and the model and reference y 2 can be written 


xV) = (y-XJ) T Y,~\y~XJ) 
Xo = y T 'P~ 1 y 


( 10 ) 

( 11 ) 


The expression for the normalized periodogram can be computed by finding via standard 
methods the value of 9 which minimizes y 2 (cu), and plugging the result into Equation ([8]). 
This yields 


Pn(u) 


F^x,. 


Tv-1 




X, 


r 1 xzx-'y 


y T ^y 


( 12 ) 


We note that this expression is equivalent to Equation ([2]) in the homoscedastic case with 
£ oc V y I. 


3.3. Simple Single-band Period Finding 


As an example of the standard periodogram in action, we perform a simple single-band 
harmonic analysis of simulated r-band observations of an RR Lyrae light curve, based on 
empirical templates derived in S10 (Figure [l]). The observations are of a star with a period 
of 0.622 days, and take place on 60 random nights over a 6-month period, as seen in the left 
panel. 


The upper-right panel shows the normalized periodogram for this source as a function 
of period. While the power does peak at the true period of 0.622 days, an aliasing effect 
is readily apparent near P = 0.38. This additional peak is due to beat frequency between 
the true period P and the observing cadence of ~ 1 day. This beat frequency is the first 
in a large sequence: for nightly observations, we’d expect to find excess power at periods 
P n = P/( 1 + nP) days, for any integer n. The strong alias in Figure [I] corresponds to the 
n — 1 beat period P n = 0.383. Though it is possible to carefully correct for such aliasing by 


iteratively removing contributions from the estimated window function (e.g. Roberts et al. 


1987), we’ll ignore this detail in the current work. 


The lower-right panel of Figure [T| shows the maximum likelihood interpretation of this 
periodogram: it is a measure of the normalized y 2 for a single-term sinusoidal model, ffere 
we visualize the data from the left panel, but folded as a function of phase, and overplotted 
with the best-fit single-term model. This visualization makes it apparent that the single¬ 
term model is highly biased: RR Lyrae light curves are, in general, much more complicated 
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Fig. 1.— An illustration of the basic periodogram and its relationship to the single-term 
sinusoid model. The left panel shows the input data, while the right panels show the fit 
derived from the data. The upper-right panel shows the periodogram with a clear peak at 
the true period of 0.622 days, and the bottom-right panel shows the data as a function of 
the phase associated with this period. Note in the periodogram the presence of the typical 
aliasing effect, with power located at beat frequencies between the true period and the 1-day 
observing cadence (see Section 3.3 for further discussion). 





than a simple sinusoid. Nevertheless, the simplistic sinusoidal model is able to recover the 
correct frequency to a high degree of accuracy (roughly related to the width of the peak) and 
significance (roughly related to the height of the peak). For a more complete introduction 
to and discussion of the single-term normalized periodogram, refer to, e.g. Bretthorst (1988) 


or 


Ivezic et al. (2014). 


4. Extending the Periodogram 

We have shown two forms of the classic normalized periodogram: Equation ([2]) and 


Equation (12). Though the two expressions are equivalent, they differ in their utility. Because 


the expression in Equation (|2]) avoids the explicit construction of a matrix, it can be computed 
very efficiently. Furthermore, through clever use of the Fast Fourier Transform, expressions 
of the form of Equation (J2| can be evaluated exactly for N frequencies in (9 [log IV] time 


(Press & Rybicki 1989). 


The matrix-based formulation of Equation (12), though slower than the Fourier-derived 


formulation, is a more general expression and allows several advantages: 


1. It is trivially extended to heteroscedastic and/or correlated measurement noise in the 
data yk through appropriate modification of the noise covariance matrix E. 

2. It is trivially extended to more sophisticated linear models by appropriately modifying 
the design matrix X^. 


3. It is trivially extended to include Tikhonov/L2-regularization terms (see Section 4.3 for 
more details) by adding an appropriate diagonal term to the normal matrix Xfl’E~ 1 X u . 


In the remainder of this section, we will explore a few of these modifications and how they 
affect the periodogram and resulting model fits. 


4.1. Stationary Sinusoid with Floating Mean 

As an example of one of these generalizations, we’ll consider what |Zechmeister & Kiirster 


(2009) call the “generalized Lomb-Scargle method”, and which we’ll call the floating-mean 
periodogram. This method adjusts the classic normalized periodogram by fitting the mean 
of the model alongside the amplitudes: 


y(t | cu, 9) = 9 0 + 9i sin cut + 9 2 cos cut 


(13) 
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Fig. 2.— An illustration of the effect of the floating mean model for censored data. The data 
consist of 80 observations drawn from a sinusoidal model. To mimic a potentially damaging 
selection effect, all observations with magnitude fainter than 16 are removed (indicated by 
the light-gray points). The standard and floating-mean periodograms are computed from 
the remaining data; these fits are shown over the data in the left panel. Because of this 
biased observing pattern, the mean of the observed data is a biased estimator of the true 
mean. The standard fixed-mean model in this case fails to recover the true period of 0.622 
days, while the floating mean model still finds the correct period. 
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This model can be more accurate than the standard periodogram for certain observing 
cadences and selection functions. Zechmeister & Kiirster (2009) detail the required modifi¬ 


cations to the harmonic formalism of Equation (J 2 j) to allow the mean to float in the model. 
In the matrix formalism, the modification is much more straightforward: all that is required 
is to add a column of ones to the X u matrix before computing the power via Equation (12). 
This column of ones corresponds to a third entry in the parameter vector 6, and acts as a 
uniform constant offset for all data points. 


For well-sampled data, there is usually very little difference between a standard peri¬ 
odogram and a floating-mean periodogram. Where this becomes important is if selection 
effects or observing cadences cause there to be preferentially more observations at certain 
phases of the light curve: a toy example demonstrating this situation is shown in Figure [2] 
The data are drawn from a sinusoid with Gaussian errors, and data with a magnitude fainter 
than 16 are removed to simulate an observational bias (left panel). Because of this observa¬ 
tional bias, the mean of the observed data are a poor predictor of the true mean, causing the 
standard fixed-mean method to poorly fit the data and miss the input period (upper-right 
panel). The floating-mean approach is able to automatically adjust for this bias, resulting 
in a periodogram which readily detects the input period of 0.622 days (lower-right panel). 


4.2. Truncated Fourier Models 


As mentioned above, the standard periodogram is equivalent to fitting a single-term 
stationary sinusoidal model to the data. A natural extension is to instead use a multiple- 
term sinusoidal model, with frequencies at integer multiples of the fundamental frequency. 
With N Fourier terms, there are 2N + 1 free parameters, and the model is given by 

N 

y(t\uj, 6) = 9 0 + [02n-i sin (nut) + 9 2n cos(naA)]. (14) 

n= 1 

Because this model remains linear in the parameters 6, it can be easily accommodated into 
the matrix formalism above. For example, an N = 2-term floating-mean model can be 
constructed by building a design matrix X u with 2N + 1 = 5 columns: 



' 1 

sin(o;fi) 

cos (cut 1 ) 

sin( 2 u+i) 

cos( 2 cnfi) 


1 

sin(o;t 2 ) 

cos ( 0 + 2 ) 

sin ( 2 ^+ 2 ) 

cos( 2 c jjt 2 ) 

V® = 

1 

sin ( 0 + 3 ) 

cos (tuts) 

sin( 2 u+ 3 ) 

cos( 2 cnf 3 ) 


_ 1 

sin (tut at) 

cos(aAjv) 

sin (2 cnt/v) 

cos( 2 a+Ar) 
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Fig. 3.— The model fits and periodograms for several truncated Fourier models. The data 
are the same as those in Figure [l] Note that the higher-order models will generally show 
periodogram peaks at multiples of the true fundamental frequency P$\ this is because for 
integer n less than the number of Fourier terms in the model, Pq is a higher harmonic of 
the model at P = nP 0 . Additionally, the increased degrees of freedom in the higher-order 
models let them fit better at any frequency, which drives up the “background” level in the 
periodogram. 
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• ( 2 ) 

Computing the power via Equation (12) using Xti ; will give the two-term periodogram. 
For larger N, more columns are added, but the periodogram can be computed in the same 
manner. Figure [3] shows a few examples of this multiterm Fourier approach as applied to 
the simulated RR Lyrae light curve from Figure [lj and illustrates several important insights 
into the subtleties of this type of multiterm fit. 


First, we see in the right panel that all three models show a clear signal at the true 
period of P 0 = 0.622 days. The higher-order models, however, also show a a spike in power 
at Pi = 2P 0 : the reason for this is that for and N > 1-term model, the period P 0 is the first 
harmonic of a model with fundamental frequency 2 P 0l and the higher-order models contain 
the single-period result. 

Second, notice that as the number of terms is increased, the general “background” level 
of the periodogram increases. This is due to the fact that the periodogram power is inversely 
related to the y 2 of the fit at each frequency. A more flexible higher-order model can better 
fit the data at all periods, not just the true period. Thus in general the observed power 
of a higher-order Fourier model will be everywhere higher than the power of a lower-order 
Fourier model. 


4.3. Regularized Models 


The previous sections raise the question: how complicated a model should we use? We 
have seen that as we add more terms to the fit, the model will come closer and closer to 
the observed data. In the extreme case, when the number of model parameters equals the 
number of data points, the model can fit the data exactly regardless of frequency and the 
periodogram will be everywhere unity (though in most cases, numerical inaccuracies prevent 
a truly perfect fit). For very high-order models, the fit becomes very sensitive to the noise 
in the data, and we say we have over-fit the data. This can be addressed by explicitly 
truncating the series, but we can also use a regularization term to mathematically enforce a 
less complicated model. 


A regularization term is an explicit penalty on the magnitude of the model parameters 
6 , and can take a number of forms. For computational simplicity here we’ll use an L2 
regularization - also known as Tikhonov Regularization (Tikhon ov|1963 ) or Ridge Regression 
(Hoerl & Kennard 11970) - which is a quadratic penalty term in the model parameters added 
to the y 2 . Mathematically, this is equivalent in the Bayesian framework to using a zero-mean 
Gaussian prior on the model parameters. 


We encode our regularization in the matrix A = diag([Ai, A 2 • • • Am]) for a model with 
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M parameters, and construct a “regularized” y 2 : 


XaM = (y - XJ) T Yr\y - XJ) + e T Ad 


(16) 


Minimizing this regularized y 2 , solving for 6, and plugging into the expression for gives 


us the regularized counterpart of Equation (12): 


Pn, a(<^) — 


y 




[XjE-^ + A]- 1 XlYT'y 


y T Z-iy 


(17) 


Notice that the effect of this regularization term is to add a diagonal penalty to the normal 
matrix X^T,~ 1 X Ul which has the additional feature that it can correct ill-posed models where 
the normal matrix is non-invertible. This feature of the regularization will become important 
for the multiband models discussed below. 


In Figure [4j we compare a regularized and unregularized 20-term truncated Fourier 
model on our simulated RR Lyrae light curve. We use A = 0 on the offset term, and 
make the penalty A j progressively larger for each harmonic component. The result of the 
regularized model is less over-fitting to the input data, and stronger periodogram peaks. 


5. A Multiple-Band Model 

In this section we will apply what was learned in the previous sections to construct 
the multiband periodogram which can flexibly account for heterogeneous sources of data 
for a single object. To compute such a periodogram, we will take advantage of the easy 
extensibility of the matrix formalism which led to our generalizations above. The multiband 
model contains the following features: 

1. An Wa-se-term truncated Fourier fit which models a latent parameter, which we’ll call 
the “overall variability”. 

2. A set of N band -term truncated Fourier fits, each of which models the residual of a single 
band from this overall variability. 


The total number of parameters for K filters is then Mk = (2 Nb ase + 1) + K(2Nf )an d + 1). 
As a result, for each band k we have the following model of the observed magnitudes: 


Vk(t\u,9) = 


/) | \ 'Mbase 

u o + Z ^„=1 

n(k) , \~^M barld 
U 0 + 1 


[d 2n -i sin (nut) + d 2n cos(no;t)] 
02 n-i sin (nut) + 6^ cos (nut) 


+ 


(18) 
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The important feature of this model is that all bands share the same base parameters 9, 
while their offsets 9^ are determined individually. 

We can construct the normalized periodogram for this model by building a sparse design 
matrix with Mk columns. Each row corresponds to a single observation through a single 
band. Columns corresponding to the base model and the matching observation band will 
have nonzero entries; all other columns will be filled with zeros. For example, the N hase = 1 
and Nband = 0 model corresponds to one with a simple single-term periodic base frequency, 
and an independent constant offset term in each band. The associated design matrix depends 
on the particular data, but will look similar to this: 



' 1 

sin(u;£i) 

cos(o;£i) 

1 

0 

0 

0 

0 ' 



1 

sin(u;f2) 

cos(o;£2 ) 

0 

0 

0 

0 
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X (l,0) = 

1 

sin(u;f3) 

cos(c at 3 ) 

0 

0 

0 

1 

0 

(19) 


_ 1 

sin(u;fjv) 

cos (tut at) 

0 

0 

1 

0 

0 . 



ffere the nonzero entries of the final five columns are binary flags indicating the (u, g, r, i, z)- 
band of the given observation: for this example, the first row is a w-band measurement, the 
second is a z-band, the third is a i-band, etc., as indicated by the position of the nonzero 
matrix element within the row. 


On examination of the above matrix, it’s clear that the columns are not linearly indepen¬ 
dent (i.e. X u is low-rank), and thus the parameters of the best-fit model will be degenerate. 
Intuitively, this is due to the fact that if we add an overall offset to the base model, this can 
be perfectly accounted for by subtracting that same offset from each residual. Mathemati¬ 
cally, the result of this is that the normal matrix will be non-invertible, and thus 

the periodogram is ill-defined. In order to proceed, then, we’ll either have to use a different 
model, or use a cleverly-constructed regularization term on one of the offending parameters. 

We’ll choose the latter here, and regularize all the band columns while leaving the base 
columns un-regularized: for the above X u matrix, this regularization will look like 

A (1,0) = diag([0, 0, 0, e, e, e, e, e]) (20) 


where e is some small fraction of the trace of the normal matrix [X^Yi~ x Xj\. The logic of 
this choice of regularization is that any component of the model which is common to all 
bands will be preferentially reflected in the base terms, with independent behavior reflected 
in the individual band terms. Setting e to some small fraction of the trace ensures that the 
regularization effect on the remaining model will be small. With this regularization in place, 
the model is well-posed and Equation (17) can be used to straightforwardly compute the 
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power. The effective number of free parameters for such a regularized ( N base , N band ) model 
with K filters is M e J f = 2N^J e + K(2N band + 1) where = max(0, N base - N band ) is the 
effective number of base terms. 

The final remaining piece to mention is our assumption in Equation (J5]) that the data 
are centered. This is required so that the simple form of the reference Xo remains valid. 
For the multiband model, this assumption requires that the data satisfy Equation (J5| within 
each band: equivalently, we could lift this assumption and compute the reference Xo °f the 
multiband model with an independent floating mean within each band; the results will be 
identical. 

This multiband approach, then, actually comprises a set of models indexed by their 
value of N base and N band . The most fundamental models have ( N base , N band ) = (1, 0) and 
(0,1), which we’ll call the shared-phase and multi-phase models respectively. In the shared- 
phase model, all variability is assumed to be shared between the bands, with only the fixed 
offset between them allowed to float. In the multi-phase model, each band has independent 
variability around a shared fixed offset. 


5.1. Relationship of Multiband and Single-band approaches 

The multi-phase (N base = 0, N band = 1) model turns out to be a particularly special case. 
Here the base model is a simple global offset which is degenerate with the offsets in each 
band, so that the design matrix X u can be straightforwardly rearranged as block-diagonal. 
A block-diagonal design matrix in a linear model indicates that components of the model 
are being solved independently: here these independent components amount to the single¬ 


band floating-mean model from Section [4Tj fit independently for each of the K bands. This 
particular multiband model can give us insight into the relationship between single-band and 
multiband approaches. 


For band k, we’ll denote the single-band floating-mean periodogram as 


AV) = i 


The full multiband periodogram is given by 


XminA^) 
Xo ,k 


( 21 ) 


E f=i xLn,fcM 




= i 


Ef=i xg,* 


( 22 ) 
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(k) 

and it can easily be shown that the Pm can be constructed as a weighted sum of : 


p(°h) 

r N 




Efc=i 

EfcLi xL 


(fc) 

JV 


(23) 


We see that this particular multiband periodogram is identical to a weighted sum of standard 
periodograms in each band, where the weights xlk are a reflection of both the number 
of measurements in each band, and how much those measurements deviate from a simple 
constant reference model. 


5.2. Multiband Periodogram for Simulated Data 


Before applying the multiband method to real data, we will here explore its effectiveness 
on a simulated RR Lyrae lightcurve. The upper panels of Figure [5] show a multiband version 
of the simulated RR Lyrae light curve from Figure [lj The upper-left panel shows 60 nights 
of observations spread over a 6-month period, and for each night all five bands ( u,g,r,i,z ) 
are recorded. Using the typical approach from the literature, we individually compute the 
standard normalized periodogram within each band: the results are shown in the upper- 
right panel. The data are well-enough sampled that a distinct period of 0.622 days can be 


recognized within each individual band, up to the aliasing effect discussed in Section 3.3 


Previous studies have made use of the information in multiple bands to choose between 
aliases and estimate uncertainties in determined periods (e.g. Oluseyi et al. [2012 Sesar et al. 


2010). While this approach is sufficient for well-sampled data, it becomes problematic when 


the multiband data are sparsely sampled. 


The lower panels of Figure [5] show the same 60 nights of data, except with only a single 
band observation recorded each night. The lower-left panel shows the observations as a 
function of phase, and the lower-right panels show the periodograms derived from the data. 
With only 12 observations for each individual band, it is clear that there is not enough data 
to accurately determine the period within each single band. The shared-phase multiband 
approach (i.e. Nb aS e = 1 , Nband = 0), shown in the lower-right panel, fits a single model to 
the full data and clearly recovers the true frequency of 0.622 days. 


This shared-phase (1,0) model is only one of the possible multiband options, however: 
Figure [6] shows multiband fits to this data for models with various choices of (N base , N band ). 
We see here many of the characteristics noted above for single-band models: as discussed in 
Section 42 increasing the number of Fourier terms leads to power at multiples of the funda¬ 
mental period, and increased model complexity (roughly indexed by the effective number of 
free parameters M e ) tends to increase the background level of the periodogram, obscuring 
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significant peaks. For this reason, models with N base > N band are the most promising: they 
allow a flexible fit with minimal model complexity. Motivated by this, in the next section 
we’ll apply the simplest of this class of models, the (1, 0) shared-phase model, to data from 
the Stripe 82 of the Sloan Digital Sky Survey. 


6. Application to Stripe 82 RR Lyrae 


Stripe 82 is a three hundred square degree equatorial region of the sky which was 
repeatedly imaged through multiple band-passes during phase II of the Sloan Digital Sky 


Survey (SDSS II, see Sesar et al. 2007). Here we consider the SDSS II observations of 483 RR 
Lyrae stars compiled and studied by S10, in which periods for these stars were determined 
based on empirically-derived light curve templates. Because the template-fitting method 
is extremely computationally intensive, S10 first determined candidate periods by taking 


the top 5 results of the Supersmoother (Reimann 1994) algorithm applied to the g- band; 


template fits were then performed at each candidate period and the period with the best 
template fit was reported as the true period. In this section, we make use of this dataset to 
quantitatively evaluate the effectiveness of the multiband periodogram approach. 


6.1. Densely-sampled Multiband Data 

The full S10 RR Lyrae dataset consists of 483 objects with an average of 55 observations 
in each of the five SDSS ugriz bands spread over just under ten years. In the upper panels 
of Figure [7] we show the observed data for one of these objects, along with the periodogram 
derived with the single-band supersmoother model and the shared-phase (0, l)-multiband 
mode0 Here we have a case which is analogous to that shown for simulated data in the top 
panels of Figure [5] each band has enough data to easily locate candidate peaks, the best of 
which is selected via the S10 template-fitting procedure. 

The lower panels of Figure [7] compare the S10 period with the best periods obtained 
from the 1-band supersmoother (lower-left) and from the shared-phase multiband model 
(lower-right). To guide the eye, the figure includes indicators of the locations of beat aliases 
(dotted lines) and multiplicative aliases (dashed lines) of the S10 period. 


lr The supersmoother “periodogram” P$s is constructed from the minimum sum of weighted model resid¬ 
uals f m i n in analogy with Equation ([8]): Pss(u) = 1 — f m i n {uj)/r o, where fo is the mean absolute residual 
around a constant model. 
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The best-fit supersmoother period matches the S10 period in 89% of cases (421/483), 
while the best-fit multiband period matches the S10 period in 79% of cases (382/483). The 
modes of failure are instructive: when the snpersmoother model misses the S10 period, it 
tends to land on a multiplicative alias (i.e. the dashed line). This is due to the flexibility 
of supersmoother: a doubled period spreads the points out, leading to fewer constraints in 
each neighborhood and thus a smaller average residual around model. In other words, the 
SuperSmoother tends to over-fit data which is sparsely-sampled. On the other hand, when 
the multiband model misses the S10 period, it tends to land on a beat alias between the S10 
period and the 1-day observing cadence (i.e. the dotted lines). This is due to the fact that 
the single-frequency periodic model is biased, and significantly under-fits the data: it cannot 
distinguish residuals due to underfitting from residuals due to window function effects. 

In both models, the S10 period appears among the top 5 periods 99% of the time: 
477/483 for supersmoother, and 480/483 for multibandj/] This suggests that had S10 used 
the multiband Lomb-Scargle rather than the supersmoother in the first pass for that study, 
the final results presented there would be for the most part unchanged. 


The results of this subsection show that the shared-phase multiband approach is com¬ 
parable to the supersmoother approach for densely-sampled multiband data, although it has 
a tendency to get fooled by structure in the survey window. Correction for this based on 
the estimated window power may alleviate this (see Roberts et al. (1987) for an example of 
such an approach) though in practice selecting from among the top 5 peaks appears to be 
sufficient. 


6.2. Sparsely-sampled Multiband Data 

Above we saw that the multiband model is comparable to methods from the literature 
for densely-sampled data. Where we expect the multiband approach to gain an advantage is 
when the data are sparsely sampled, with data through only a single band at each observation 
time. To simulate this, we reduce the size of the Stripe 82 RR Lyrae dataset by a factor of 
5, keeping only a single band of imaging each night: an average of 11 observations of each 
object per band. This is much closer to the type of data which will be available in future 
multiband time-domain surveys. 


2 We might expect this correspondence to be 100% in the case of the g-band supersmoother, which was 
the model used in the first pass of the S10 computation. This discrepancy here is likely due to the slightly 
different supersmoother implementations used in S10 and in this work. Objects showing this discrepancy 
are those with very low signal-to-noise. 
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The upper panels of Figure [8] show an example light curve from this reduced dataset, 
along with the supersmoother and multiband periodograms derived from this data. Analo¬ 
gously to the lower panels of Figure [5j the single-band supersmoother model loses the true 
period within the noise, while the shared-phase multiband model still shows prominent signal 
near the S10 period. 

The lower panels of Figure [8] show the relationship between the S10 periods (based on 
the full dataset) and the periods derived with each model from this reduced dataset. It is 
clear that the supersmoother model is simply over-fitting noise with this few data points: 
the top period matches S10 in only 32% of cases (compared to 87% with the full dataset), 
and the top 5 periods contain the S10 period only 45% of the time. The failure mode is 
much less predictable as well: rather than being clustered near aliases, most of the period 
determinations are scattered seemingly randomly around the parameter space. 

The multiband method does much better on the reduced dataset. Even with an 80% 
reduction in the number of observations, the multiband method matches the S10 period 
64% of the time (compared to 79% with the full dataset), and the top 5 peaks contain the 
S10 period 94% of the time (compared to 99% with the full dataset). This performance is 
due to the fact that the multiband algorithm has relatively few parameters, but is yet able 
to flexibly accommodate data from multiple observing bands. In particular, this suggests 
that with the multiterm periodogram, the S10 analysis could have been done effectively with 
only a small fraction of the available data. This bodes well for future surveys, where data 
on variable stars will be much more sparsely sampled. 


6.3. Potential Improvements to the Multiband Method 


The fact that the multiband model does not select the top frequency each time points 
to its weaknesses: first, as a Fourier-like model, it tends to respond to frequency structure in 
the window function as well as frequency structure in the data. This is a result of the very 
model simplicity which causes its success in the case of sparse multiband data: it cannot 
disentangle bias in the model from bias due to the survey window. This could potentially 
be accounted for by correcting for the effect of the estimated window function; one potential 
method for this involves estimating the deconvolution of the window power and the observed 


power (Roberts et al. 1987). It may also be possible to propose a multiband extension of 


CARMA (Kcllyet al. 2014), or another forward-modeling approach to detecting periodicity. 


One potentially fruitful avenue of research which we do not study here is the application 
of other types of regularization to the higher-order periodogram. In particular, LI regular- 
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ization (also known as Lasso regression) could lead to interesting results: LI regularization 


is similar in spirit to the Tikhonov regularization discussed in Section 4.3, but tends toward 


sparsity in the model parameters (see, e.g. Ivezic et al. 2014) for a discussion). Such an 
approach could provide a useful tradeoff between model complexity and bias in the case of 
higher-order truncated Fourier models. 

Another potentially interesting extension of the multiband case would be to define and 
make use of physically-motivated priors in the light-curve shape. This approach could allow 
the model bias to be decreased without a commensurate increase in model complexity, which 
is what causes poor performance in the case of sparsely-sampled noisy data. As an example of 
such a physically-motivated prior, consider that the paths of RR Lyrae stars through color- 
color and color-magnitude space are constrained by known astrophysical processes in the 


structure of the stars (e.g., see Fig. 5 in Szabo et al. 2014). Making use of this information 
could help break degeneracies in period determination with higher-order models. 


7. Prospects for Multiband Periodograms with LSST 


Previously, Oluseyi et al. (2012) evaluated the prospects of period finding in early LSST 
data, and found results which were not encouraging. Using the conservative criterion of a 
2/3 majority among the top single-band supersmoother periods in the g, r, and i bands, they 
showed that, depending on spectral type, finding reliable periods for the brightest (g ~ 20) 
RR Lyrae stars will require several years of LSST data, while periods for some of the faintest 
(g ~ 25) stars will not be reliable with even ten years of data! 


One potential remedy is to move away from general models like supersmoother and 
lomb-scargle to specific template-fitting methods such as those used in S10. Indeed, such 
methods perform well even for sparsely-sampled multiband data such as those from the 
PanSTARRS survey; the primary drawback is that such blind template fits are computation¬ 
ally extremely expensive: they involve nonlinear optimizations over each of several hundred 
candidate templates at each of tens of thousands of candidate frequencies (B. Sesar, private 
communication). Thus the template-fitting method, though it produces accurate periods, 
in practice requires several hours of CPU time for a well-sampled period grid for a single 
source (compared to several seconds for the multiband periodogram proposed here). Note 
that several hours per object is orders-of-magnitude too slow in the case of LSST; to estimate 
periods for a billion stars on a 1000-core machine in a year requires a compute-time budget 
of only 30 seconds per light curve. 


Because of the computational expense of the pure template-fitting method, when work- 
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ing with SDSS II data S10 performed a first-pass with a single-band supersmoother to 
establish candidate periods, which were in turn evaluated with template-fitting approach. 
Here we show that such a hybrid strategy combining the multiband periodogram and the 
S10 template fits will be useful for determining periodicity of variables in early LSST data 
releases, greatly improving on the outlook presented in Oluseyi et ah (2012). 

We suggest the following procedure for determining periods in future multiband datasets: 


1. As a first pass, find a set of candidate frequencies using the multiband periodogram. 
This is a fast linear optimization that can be trivially parallelized. 

2. Within these candidate frequencies, use the more costly template-fitting procedure to 
choose the optimal period from among the handful of candidates. 

3. Compute a goodness-of-fit statistic for the best-fit template to determine whether the 
fit is suitable; if not, then apply the template-fitting procedure across the full period 
range. 

Here we briefly explore simulated LSST observations of RR Lyrae stars in order to gauge 
the effectiveness of the first step in this strategy; the effectiveness of the template-fitting 
step will be explored further in future work. Rather than doing the full analysis including 
the final template fits, we will focus on the ability of the multiband periodogram to quickly 
provide suitable candidate periods under the assumption that the S10 template algorithm 
will then select or reject the optimal period from this set. 


7.1. LSST Simulations 


We use a simulated LSST cadence (Delgado et al.||2006 Ridgway et al.||2012 Jones et al. 


2014) in 25 arbitrarily chosen fields that are representative of the anticipated main survey 


temporal coverage. We simulate a set of 50 RR Lyrae observations with the S10 templates, 
with a range of apparent magnitudes between g = 20 and g = 24.5, corresponding to bright- 
to-faint range of LSST main-survey observations, and with expected photometric errors 


computed using eqs. 4-6 from Ivezic et al. (2008). Given the capability of template-fitting 


to choose among candidate periods, we use a more relaxed period-matching criterion than 


m 


Oluseyi et al. (2012): when evaluating the single-band supersmoother, we require that 


the true period is among the five periods determined independently in the u , g , r, i , z bands; 
in the multiband case we require that the true period is among the top five peaks in the 
multiband periodogram. 
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Figure [9] shows the fraction of stars where this period matching criterion is met as a 
function of g-band magnitude and subset of LSST data. The solid lines show the multiband 
results; the dashed lines show the single-band supersmoother results; and the shading helps 
guide the eye for the sake of comparison. Because of our relaxed matching criteria, even 


the single-band supersmoother results here are much more optimistic than the Oluseyi et al. 


(2012) results (compare to Figure 15 in that work): the supersmoother result here can be 
considered representative of a best-case scenario for ad hoc single-band fits. Without fail, the 
multiband result exceeds this best-case single-band result; the improvement is most apparent 
for faint stars, where the greater model flexibility of the supersmoother causes it to over-fit 
the noisy data. 


The performance of the multiband periodogram points to much more promising prospects 
for science with variable stars than previously reported. In particular, even with only six 
months of LSST data, we can expect to correctly identify the periods for over 60% of stars 
brighter than g = 22; with the first two years of LSST observations, this increases to nearly 
100%; with five years of data, the multiband method identifies the correct period for 100% 
of even the faintest stars. Part of this improvement is due to the performance of the shared- 
phase multiband model with noisy data, and part of this improvement is due to the re¬ 
laxed period-matching constraints enabled by the hybrid approach of periodogram-based 
and template-based period determination. 


8. Discussion and Conclusion 

We have motivated and derived a multiband version of the classic Lomb-Scargle method 
for detecting periodicity in astronomical time-series. Experiments on several hundred RR 
Lyrae stars from the SDSS Stripe 82 dataset indicate that this method outperforms meth¬ 
ods used previously in the literature, especially for sparsely-sampled light curves with only 
single bands observed each night. While there are potential areas of improvement involv¬ 
ing corrections to window function artifacts and accounting for physically-motivated priors, 
the straightforward multiband model outperforms previous ad hoc approaches to multiband 
data. 

Looking forward to future variable star catalogs from PanSTARRS, DES, and LSST, 
there are two important constraints that any analysis method must meet: the methods must 
be able to cope with heterogeneous and noisy observations through multiple band-passes, 
and the methods must be fast enough to be computable on millions or even billions of 
objects. The multiband method, through its combination of flexibility and model simplicity, 
meets the first constraint: as shown above, in the case of sparsely-sampled noisy multiband 
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data, it out-performs previous approaches to period determination. It also meets the second 
constraint: it requires the solution of a simple linear model at each frequency, compared to 
a rank-based sliding-window model in the case of supersmoother, a nonlinear optimization 
in the case of template-fitting, and a Markov Chain Monte Carlo analysis in the case of 
CARMA models. In our own benchmarks, we found the multiband method to be several 
times faster than the single-band supersmoother approach, and several orders of magnitude 
faster than the template fitting approach. 

The strengths and weaknesses of the multiband method suggest a hybrid approach to 
finding periodicity in sparsely-sampled multiband data: a first pass with the fast multiband 
method, followed by a second pass using the more computationally intensive template-fitting 
method to select among these candidate periods. Despite pessimism in previous studies, our 
experiments with simulated LSST data indicate that such a hybrid approach will successfully 
identify periods in the majority of RR Lyrae stars brighter than g ~ 22.5 in the first months 
of the survey, and the majority of the faintest detected stars with several years of data. This 
finding suggests that the multiband periodogram could have an important role to play in 
the analysis of variable stars in future multiband surveys. 


We have released a Python implementation of the multiband periodogram on GitHub, 
along with Python code to reproduce all results and figures in this work; this is described in 
Appendix [Aj As we were finalizing this manuscript, we were made aware of a preprint of an 
independent exploration of a similar approach to multiband light curves (Long et ah 2014); 


we 


discuss the similarities and differences between these two approaches in Appendix |B| 
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A. Python Implementation of Multiband Periodogram 


The algorithm outlined in this paper is available in gat spy, an open-source Python 
package for general astronomical time-series analysi^] (Vanderplas 2015a). Along with the 
periodogram implementation, it also contains code to download all the data used in this work. 
Code to reproduce this paper, including all figures, is available in a separate repositor)]^} 

gat spy is a pure-Python package written to be compatible with both Python 2 and 


Python 3, and performs fast numerical computation through dependencies on numpy (van der 


Walt et al. 

2011 

3 and astroML ( 

Vanderplas et al. 

2012 6 


which offer optimized implemen¬ 


tations of numerical methods in Python. 


The API for the module is largely influenced by that of the scikit-learn package 

in which models are Python class objects 


(Pedregosa et al. 

2011 

Buitinck et al. 

2013 7 


which can be fit to data with the f it () method. Here is a basic example of how you can use 
multiband_LS to download the data used in this paper, fit a multiband model to the data, 
and compute the power at a few periods: 


from gatspy.periodic import LombScargleMultiband 
import numpy as np 


# Fetch, the Sesar 2010 RR Lyrae data 
from gatspy.datasets import fetch.rrlyrae 
data = fetch.rrlyrae() 

t, mag, dmag, filts = data.get_lightcurve(data.ids [0]) 


# Construct the multiband model 

model = LombScargleMultiband(Nterms_base=0, Nterms_band=l) 
model.fit(t, mag, dmag, filts) 


# Compute power at the following periods 

periods = np.linspace(0.2, 1.4, 1000) # periods in days 

power = model.periodogram(periods) 


Other models are available as well. For example, here is how you can compute the 


a http://github.com/astroml/gatspy/ 

J http://github.com/jakevdp/multiband_LS/ 


http: //www. numpy. org 


<: http: //www. astroml. org 

'http://scikit-learn. 

org 
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periodogram under the supersmoother model; this implementation of the supersmoother 
periodogram makes use of the supersmoother Python package (Vanderplas 2015b). 

from gatspy.periodic import SuperSmoother 

# Construct the supersmoother model 
model = SuperSmoother() 

gband = (filts == J g’) 

model.fit(t[gband] , mag[gband] , dmag[gband]) 

# Compute power at the given periods 
power = model.periodogram(periods) 


(Vanderplas 2015b 


The models in the gatspy package contain many more methods, and much more func¬ 
tionality that what is shown here. For updates, more examples, and more information, visit 
http://github.com/astroml/gatspy/. 


B. 


Comparison with Long et al. (2014) 


As we were finishing this study, we learned that another group had released a preprint 
independently addressing the multiband periodogram case, and come up with a solution very 
similar to the one presented here (Long et al. 2014, hereafter LCB14). They present two 
methods, the “Multiband Generalized Lomb-Scargle” (MGLS) which is effectively identical 
to the (1, 0) multi-phase model here, and the “Penalized Generalized Lomb-Scargle” (PGLS), 
which is similar in spirit to our (0,1) shared-phase model. 


In the PGLS model, they start with a multi-phase model, fitting independent N = 1 
term fits to each band, and apply a nonlinear regularization term which penalizes differences 
in the amplitude and phase. In terms of the formalism used in this work, the PGLS model 
minimizes a regularized y 2 of the form 


xIcls = E [ xIls(D { k) ) + Ja(AM) + JM {k) ) 

k =1 


(Bl) 


for K bands, where Xgls( is the x 2 °f the standard floating mean model on the single¬ 
band data ZA A ), and Ja and J$ are regularization/penalty terms which are a function of the 
amplitude A k and phase (p^ k> of each model. In terms of our linear model parameters 0 ( ~ k \ 
this amplitude and phase can be expressed: 

= \J (O^) 2 + (# 2 fe) ) 2 

<j>^ — arctan($ 2 /$i) 


(B2) 












The selected form of these regularization terms penalizes deviations of the amplitude and 
phase from a common mean between the bands; in this sense the PGLS model can be consid¬ 
ered a conceptual mid-point between our shared-phase and multi-phase models. Within the 
formalism proposed in the current work, such a mid-point may be alternatively attained by 
suitably increasing the regularization parameter e used in our shared-phase model, though 
the nature of the resulting regularization will differ. 

Computationally, the PGLS model requires a nonlinear optimization at each frequency 
oj, and is thus much more expensive than the straightforward linear optimization of our 
shared-phase model. For this reason, LCB14 proposes a clever method by which nested 
models are used to reduce the number of nonlinear optimizations used: essentially, by showing 
that the (linear) MGLS y 2 is a lower-bound of the (non-linear) PGLS y 2 , it is possible to 
iteratively reduce the number of PGLS computations required to minimize the y 2 among a 
grid of frequencies. Such an optimization could also be applied in the case of our shared-phase 
model, but is not necessary here due to its already high speed. Nevertheless, when applying 
the method to a very large number of light curves, as in e.g. LSST, such a computational 
trick may prove very useful. 

Given these important distinctions between the models proposed here and in LCB14, 
in future work we plan to do a detailed comparison of the two means of multiband model 
regularization. 
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Fig. 4.— The effect of regularization on a high-order model. The data is the same as 
those in Figure [lj We fit a 20-term truncated Fourier model to the data, with and without 
a regularization term. Without regularization, the model oscillates widely to fit the noise 
in the data. The regularization term effectively damps the higher-order Fourier modes and 
removes this oscillating behavior, leading to a more robust model with stronger periodogram 
peaks. 
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Fig. 5.— An illustration of the performance of the multiband periodogram. The upper 
panels show simulated ugriz observations of an RR Lyrae light curve in which all 5 bands 
are observed each night. With 60 observations in each band, a periodogram computed from 
any single band is sufficient to determine the true period of 0.622 days. The lower panels 
show the same data, except with only a single ugriz band observed each night (i.e. 12 
observations per band). In this case, no single band has enough information to detect the 
period. The shared-phase multiband approach of Section [5] (lower-right panel) combines the 
information from all five bands, and results in a significant detection of the true period. 
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Periodograms for Multiterm Models 
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Fig. 6.— Comparison of the periodograms produced by various multiband models. The 
data is the same as that used in Figure [5j N baS e gives the number of Fourier terms in the 
base model, and N band gives the number of Fourier terms used to fit the residuals around 
this model within each band. The characteristics discussed with previous figures are also 
seen here: in particular, the level of “background noise” in the periodogram grows with the 
model complexity M, 
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Fig. 7.— Comparison of the Multiband algorithm and single-band supersmoother algorithm 
on 483 well-sampled RR Lyrae light curves from Stripe 82. The upper panels show a rep¬ 
resentative lightcurve and periodogram fits, while the bottom panels compare the derived 
periods to the template-based periods reported in S10. Shown for reference are the beat 
aliases (dotted lines) and the multiplicative alias (dashed lines): numbers along the top and 
right edges of the panels indicate the number of points aligned with each trend. The single¬ 
band supersmoother model tends to err toward multiplicative aliases, while the multiband 
model tends to err toward beat frequency aliases. Both methods find the correct period 
among the top 5 significant peaks around 99% of the time. 
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Fig. 8.— This figure repeats the experiment shown in Figure [7] (see caption there for 
description), but the data is artificially reduced to only a single-band observation on each 
evening, a situation reflective of the observing cadence of future large-scale surveys. In this 
case, the single-band SuperSmoother strategy used as a first pass in S10 fails: there is simply 
not enough data in each band to recover an accurate period estimate. The correct period 
is among the top 5 candidates in fewer than 50% of cases. The shared-phase multiband 
approach utilizes information from all five bands, and returns much more robust results: 
even with the greatly-reduced data, the true period is among the top 5 candidates in 93% 
of cases. 
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Fig. 9.— Fraction of periods correctly determined for LSST RR Lyrae as a function of the 
length of the observing season and the mean g-band magnitude, for the multiband (solid 
lines) and single-band supersmoother (dashed lines) approaches. The multiband method 
is superior to the single-band supersmoother approach in all cases, and especially for the 
faintest objects. 





